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The dynamics of a many-body system coupled to an external environment represents a funda- 
mentally important problem. To this class of open quantum systems pertains the study of en- 
ergy transport and dissipation, dephasing, quantum measurement and quantum information theory, 
phase transitions driven by dissipative effects, etc. Here, we discuss in detail an extension of time- 
dependent current-density-functional theory (TDCDFT), we named stochastic TDCDFT [Phys. 
Rev. Lett. 98, 226403 (2007)], that allows the description of such problems from a microscopic 
point of view. We discuss the assumptions of the theory, its relation to a density matrix formalism, 
and the limitations of the latter in the present context. In addition, we describe a numerically con- 
venient way to solve the corresponding equations of motion, and apply this theory to the dynamics 
of a ID gas of excited bosons confined in a harmonic potential and in contact with an external bath. 
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I. INTRODUCTION 



Density functional theory (DFT)^'^ 



has found 



widespread application in different fields ranging from 
materials science to biophysics. Its original formula- 
tion dealt with the ground-state properties of many- 
particle systems, but since then it has been extended 
to the time domain;^!^ giving access to relevant infor- 
mation about the non-equilibrium properties of many- 
body systemsi^ According to which variable is employed 
as the basic physical quantity of interest, namely the den- 
sity or the current density, these dynamical extensions 
are named time-dependent DPT (TDDFT)^ or time- 
dependent current-DFT (TDCDFT) The successes 
of these theories are impressive and are mainly due to 
their conceptual and practical simplicity which allows the 
mapping of the original interacting many-body problem 
into an effective single-particle problem. From a compu- 
tational point of view this represents a major simplifica- 
tion compared to other, equally valid, but computation- 
ally more demanding many-body techniques. 

Nevertheless, one needs to recognize that in its present 
form DFT can only deal with systems evolving under 
Hamiltonian dynamics. This leaves out a large class of 
physical problems related to the interaction of a quan- 
tum system with one or several external environments, 
namely the study of the dynamics of open quantum 
systems. ^'^'^ Examples of such problems include energy 
transport driven by a bath (e.g., thermoelectric effects), 
decoherence, phase transitions driven by dissipative ef- 
fects, quantum information and quantum measurement 
theory, etc. The study of these problems from a micro- 
scopic point of view would give unprecedented insight 
into the dynamics of open quantum systems. 

The present authors have recently extended DFT 
to the study of the dynamics of open quantum sys- 
tems by proving that, given an initial condition and a 
set of operators that describe the system-bath interac- 
tion, there is a one-to-one correspondence between the 
ensemble-averaged current density and the external vec- 



tor potential.^^ This theory has been named stochastic 
time-dependent current-DFT (S-TDCDFT)^ Its start- 
ing point is a stochastic Schrodinger equation (SSE)i 
which describes the time-evolution of the state vector in 
the presence of a set of baths, which introduce stochastic- 
ity in the system dynamics at the Markov-approximation 
level, or if the baths' operators depend locally on time, it 
represents a form of non-Mar kovian dynamics, whereby 
the interaction of the baths with the system changes in 
time, but it carries information only at the time at which 
the state vector is evaluated, and not on its past dynam- 
ics [see Eq. (g])]. A practical application of S-TDCDFT to 
the decay of excited He and its connection with quantum 
measurement theory can be found in Ref. 11. 

If the Hamiltonian of the system does not depend on 
microscopic degrees of freedom, such as the density or 
the current density, the SSE is the stochastic unraveling 
of a quantum master equation for the density matrixi^ii^ 
One could thus argue that an equation of motion for the 
many-body density matrix is an equally valid starting 
point for a functional theory of open quantum systems 
Unfortunately, this is not the case for several reasons. 
These are mainly related to the lack of a closed equation 
of motion for the density matrix when the Hamiltonian 
of the system depends on microscopic degrees of free- 
dom, and the possible lack of positivity of the density 
matrix when the Hamiltonian and/or bath operators are 
time-dependent: the Kohn-Sham (KS) Hamiltonian is, 
by construction, always time-dependent in TDDFT. As 
we will discuss in this paper, these fundamental draw- 
backs do not pertain to the solution of the SSE, making 
it a solid starting point to develop a stochastic version of 
DFT. 

The paper is organized as follows. In Sec. |TT] we in- 
troduce the basic notation of stochastic processes and 
equations of motion. In Sec. IIIII we discuss S-TDCDFT 
and in Sec. IIVI we make a connection with a density- 
matrix approach, showing the limitations of the latter in 
the present DFT context. In Scc.|V]we describe numeri- 
cally convenient ways to solve the equations of motion of 
S-TDCDFT, and in Sec. EJ we apply this theory to the 
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time evolution of a gas of excited bosons confined in a 
harmonic potential, interacting at a mean-field level and 
coupled to an external time- independent environment. 
We finally report our conclusions and plans for future 
directions in Sec. IVIII 



II. BASIC NOTATION 

Let us consider a quantum-mechanical system of N 
interacting particles of charge e subject to an external 
deterministic perturbation. The Hamiltonian of this sys- 
tem is 



H 



N 
1=1 



2m 



N 

+ Y.U^nt{h~fJ) (1) 



where Af,xt {f, t) is the external vector potential and 
Uintir) describes the particle-particle interaction poten- 
tial. We work here in a gauge in which the scalar poten- 
tial is set to vanish identically. 

Let us assume that this quantum-mechanical system is 
coupled, via given many-body operators, to one or many 
external environments that can exchange energy and mo- 
mentum with the system. If we assume that the dynam- 
ics of each environment is described by a series of inde- 
pendent memory-less processes, the dynamics of the sys- 
tem is governed by the stochastic Schrodinger equation^ 
(h = I throughout the paper) 



a 

+ j2ut)Va{mit)) 



(2) 



where [/„ and Va describe the coupling of the system 
with the a-th environment. We will see below that, if we 
impose that the state vector has an ensemble-averaged 
norm equal to one, then Ua = V^Va [Eq. IHl)], which 
provides an intuitive interpretation of these two operators 
in terms of dissipation and fluctuations, respectively. 

One can postulate that such stochastic equation gov- 
erns the dynamics of our open quantum system^i^ or, if 
the Hamiltonian is noi stochastic (i.e., it does not depend 
on microscopic degrees of freedom such as the density or 
current density), the SSE Q can be justified a posteri- 
ori by proving that it gives the correct time evolution of 
the many-particle density matrix, namely it is the un- 
raveling of a quantum master equation for the density 
matrix (see also Sec. IIVI) .^ Or better yet, one can derive 
the SSE ^ from first principles using, e.g., the Fesh- 
bach projection-operator method to trace out (from the 
total Hamiltonian: system plus environment (s) and their 
mutual interaction) the degrees of freedom of the envi- 
ronment (s) with the assumption that the energy levels of 
the latter form a dense setfi^ In this way, one can in fact 
derive an equation of motion more general than the SSE 
(12) which is valid also for environments that do not fulfill 



the memory-less approximation. In the memory-less ap- 
proximation that equation of motion reduces to the SSE 

Here we do not restrict the theory to time-independent 
Ua and Va operators, but we assume that the dynam- 
ics of these operators is not affected by the presence of 
the quantum mechanical system, i.e., we neglect possible 
feedback of the quantum mechanical system on the ex- 
ternal environments. Moreover, we assume that Va and 
Ua admit a power expansion in time at any timei^ For 
instance, a sudden switch of the system-bath coupling 
cannot be treated in our formalism. Finally we admit 
that Ua and Va may vary in space. In the following the 
time and spatial arguments of Ua and Va are suppressed 
to simplify the notation. 

We choose {Ua} to be Hermitian operators. Indeed, 
any anti-hermitian part of the Ua operators is effectively 
an external non-dissipative potential that can be included 
in the Hamiltonian, and then via a gauge transformation 
in the vector potential. In Eq. ([2]), {la{t)} are a set of 
Markovian stochastic processes 



LWpW) = Sa,l36{t - t') 



(3) 
(4) 



where the symbol indicates the stochastic average over 
an ensemble of identical systems evolving according to 
the stochastic Schrodinger equation ([2]). 



A. Ito calculus 

Clearly, Eq. ^ does not follow the "standard" rules 
of calculus. Indeed, since \^{t)) is a stochastic function 
of time its time derivative is not defined at any instant 
of time, namely, the stochastic terms, ZQ,(i)V^ and the 
Markov approximation, Eq. @ , make this equation non- 
tractable with the standard calculus techniques.* In par- 
ticular, one has to assign a meaning to quantities like 



f{t')la{t')dt' = f f{t')dWa{t') (5) 



where f{t) is a test function and Wa{t) is a Wiener pro- 
cess such thati 



Wait) = / la{t')dt 



(6) 



There are many different ways to assign a physical and 
mathematical interpretation to Eq. ([5]). In this paper we 
use the Ito calculus^ 

/ f{t')dWa{t') ^ lim y^f{t,)[Wa{t^+l)-Wa{td] 

Jo Q^°° ~{ 

(7) 

where {ti} is a series of time steps such that ti — Q 
and tn = t. For instance, another possible choice is 
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(Stratonovich) 



that is as an infinitesimal difference equation 
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/ f{t)dW^{t)^ \im 

Jo Q^oo 



fiU) + fiU+i) 



2 (8) 

x[Wo.{U+i)-W^iU)]. 

In standard calculus, one can prove that the r.h.s. of 
Eqs. ^ and ([8]) are identical. However, this is not true 
if Wa describes a stochastic process: Eqs. (O and ^ 
bear different physical interpretations, and it is then not 
surprising that they do not coincide. 

The Wiener process Wa describes the dynamics of the 
fluctuations due to the environment and defines the cou- 
pling between these fiuctuations and the system. In con- 
sidering the cumulative effect of these fiuctuations on the 
system we have (at least) two possible choices. On the 
one hand, we may assume that the only knowledge [em- 
bodied by the function f{t) in ([7]) and ([8])] on the system 
we have access to is that at times preceding the instant at 
which a fluctuation takes place, thus leading to Eq. ([7]). 
Alternatively, we can assume that the response of the 
system is determined by its properties "in between" the 
states before and after the fluctuation has occurred, and 
thus Eq. ([5]) follows. This second interpretation is correct 
only if the fluctuations of the environment are "regular" , 
i.e., if the r.h.s. of the Eq. (|4]) is replaced by a regular 
function oi t ~ t' . We will however, restrict ourselves to 
the case in which Eq. Q is valid. This has some mathe- 
matical advantages, and it is always possible to transform 
the results from one formalism to the other by a simple 
mapping.— 



B. Stochastic Schrodinger equation 

Once we have defined the rules of integration with re- 
spect to the Wiener process, the SSE ^ has to be inter- 
preted as 



l*>, (9) 



It is important to bear in mind that if the Ito approach 
is used, few of the rules of the standard calculus have to 
be modified. The most important and relevant for our 
following discussion is the rule of product differentiation 
or chain ruleM^ Indeed, we have that if ^ and are 
two states evolving according to the SSE then 



(10) 



When Eq. ^ is used to express Eq. (fTO|) in terms of the 
Hamiltonian, the following simple rules of calculus must 
be kept in mind^^ 



dtdt = 0; dtdWa = 0; dW^dW/s = 5a,(}dt. (11) 
These relations, that we assume here valid without fur- 
ther discussion, can be proved exactly in the Ito ap- 
proach to stochastic calculus.— The first two mean that 
terms of order higher than dt are neglected [from Eqs. @ 
and (III) we see that dWa ~ Vdi] while the third ensures 
that the different environments act independently on the 
dynamics of the quantum-mechanical system. 

Eqs. (fTO|) and (fTT|) will be used as basic rules of cal- 
culus throughout this paper. To simplify the notation, 
in the following we will consider only one environment. 
The generalization to many independent environments is 
straightforward. 

Having set the mathematical rules, we can now de- 
rive the equations of motion for the particle density and 
current density. These equations of motion will be our 
starting point to develop stochastic TDCDFT. By using 
Ito formula Eq. (|10p we immediately obtain the equation 
of motion for the many-particle density (this is a function 
of N coordinates, including spin) 



J 



d(***) = i(**iJ)* - i^*{H'^) - + (**l/t)(i/5,) 

I 



dt 



dW. 



(12) 



By integrating over all degrees of freedom of all particles, 
and taking the ensemble average of the result, we obtain 
the equation of motion for the ensemble-averaged total 
norm, N 



dN 
It 



{V^V-^U), 



(13) 



where the symbol (A) indicates the standard quantum- 



mechanical expectation value of the operator A. From 
Eq. (fT3|l we immediately see that if we assume = U 
we obtain that the state vector has an ensemble-averaged 
constant norm. In the following, we are then going to 
assume that 

V^V = U. (14) 
This relation is reminiscent of the "fluctuations- 
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dissipation theorem" which relates the dissipation that 
drives the system towards an equihbrium state [the terms 
i (Ja dt in Eq. with the fluctuations induced by 
the external environment (the terms VadWa in the 
same equation) and which drive the system out of equi- 
librium. Here, however, this relation is not limited to a 
system close to equilibrium but it pertains also to sys- 
tems far from equilibrium. 

Using Eq. ([H]), Eqs. ® and ^ simphfy to (for one 
environment) 



d\^) = 



dt + V\^)dW, (15) 



and 



d(^'**) = i(«'*i/)^'-i1'*(i?*)-^'*y^l/* 
+ {^*V^)V^] dt 



(16) 



(**!>)«' + **(!/*) dW, 



respectively. Starting from Eqs. (fT5)) and we can 
obtain the equation of motion for the expectation value 
of any observable A 

d{A) = (d(*|)i|*) + (*|i(d|*)) + (d(*|)i(d|*)) 

i[A, H]-^ (v'lVA + AV'<V - 2V^AV^ \ dt 



+ {V^A + Av)dW. 



(17) 



The equation of motion for the ensemble-averaged expec- 
tation value is obtained immediately from (jl7p . 



dt{A) = ^{[A,H]) 



-- { (VVA) + (AViV) 



A,H 



2{V'fAV)) (18) 



((ytf i) + (AV^'V) - 2(T>Ul>)) ,(19) 

where we have used dW = 0. 

In the last step we have also assumed that {[A,H]) = 

{[A,H]). This relation is valid only if H does not de- 
pend on any stochastic field, i.e., it is not a stochastic 
Hamiltonian which is different for the different elements 
of the statistical ensemble (see Fig. [1]). If, for example, 
the particle-particle interaction in H is treated in the 
Hartree approximation, then the last step in (jl9p is not 
justified, and the equation of motion for the expectation 
value of any operator A will not be given by Eq. ([19]) but 
by the more complex Eq. (US]). 



C. Quantum master equation 

For the simpler case in which the Hamiltonian is not 
stochastic one can easily obtain a closed equation of mo- 
tion for the density matrix from the SSE. Quite generally 



we define 

m^JwWWi ^ E (20) 

i 

where |5'i(i)) is a pure state vector in the Hilbert space 
of the system occurring in the ensemble with probability 
Pi{t), with J2i Pi{t) = 1- Definition (|20p is valid when 
the initial state of the system is pure. If the initial state 
of the system is mixed with macro-state {|^'o)iPn}i then 
definition (|20p of statistical operator must include an ex- 
tra summation 



p(t)-E^'«l*"W)(*"WI' 



(21) 



where |*"(i)) = {\^f{t))} is the ensemble of state vec- 
tors corresponding to the initial condition l^'g). Equa- 
tion (PT|) reduces to (^0)1 for a pure initial state {p^} = 
{1,0,. ..,0}. 

Using the definition (pij) of density matrix we can de- 
fine the ensemble average of any observable A as 



(i) = Tr{/5(<)i}. 



(22) 

By using Eq. (fT9]) which is valid for any observable, the 
many-particle density matrix operator follows the equa- 
tion of motion 



dtp = 



P, 



H 



pV^V -2VpvA (23) 



which is the well-known quantum master equation (or 
Lindblad equation if all operators, including the Hamil- 
tonian, do not depend on time)^^!^^ 

We stress once more that, in order to derive this quan- 
tum master equation, we have assumed that the Hamil- 
tonian does not depend on any stochastic field. Other- 
wise, our starting point would have been Eq. (fTS)) and no 
closed equation of motion for the density matrix could 
have been obtained. 

Note that this is true even if the system does not 
interact with an external environment but its state is 
mixed. A stochastic Hamiltonian prevents us from writ- 
ing a closed equation of motion for the density matrix 
while the SSE (fT^ contains this case quite naturally: one 
simply evolves the system dynamics over the ensemble of 
stochastic Hamiltonians and then averages the resulting 
dynamics. This point is particularly relevant in DFT 
where the KS Hamiltonian does depend on microscopic 
degrees of freedom, and it is thus generally stochasticii^ 

There is another important reason for not using the 
quantum master equation (|23p in a DFT approach. In 
fact, it is only when the Hamiltonian of the system and 
the bath operators are time-independent that one can 
prove that the density matrix solution of Eq. (|^5)) fulfills 
the usual requirements of a "good" statistical operator, 
i.e., that at any instance of time its trace is conserved, 
the operator is hermitian and that it remains a definite- 
positive operator, namely for any state $ in the Hilbert 
space 



($1/3(01$) > 0. 



(24) 
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operators A and B. From Eq. p8l) we then get 



dn{r, t) 



dt 



+ ^(2V^h{r, t)V - yt\/fi(r, t) ~ n{r,t)V'fV 



(29) 



FIG. 1; (Color online) If the Hamiltonian depends on mi- 
croscopic degrees of freedom such as the particle or current 
density, it is different for each element of the statistical en- 
semble represented by the probabilities pi that the system 
is found in the state vector of M accessible states. A 
stochastic Hamiltonian precludes the derivation of a closed 
equation of motion for the density matrix [see discussion fol- 
lowing Eq. (fT9)l ]. 



The reason for these restrictions is because a dynamical 
semi-group (in the exact mathematical sense) can only 
be defined for time-independent HamiltoniansJ^iiiiiSiii^ 

It is important to realize that an approach based on the 
SSE (fT5|) does not suffer from this drawback: the density 
matrix (j2ip constructed from the SSE is by definition 
positive at any time. 

All of this points once more to the fact that Eq. ([^5)1 
is not a good starting point to build a stochastic version 
of DFT. We will expand a bit more on these issues in 
Sec. IIVI In Sec. IVII we will provide an explicit example 
that shows that Eq. leads to the wrong dynamics in 
the presence of interactions among particles. 



D. Continuity equation 

We can use the general result Eq. to derive the 
equation of motion for the ensemble-averaged particle 
density. Let us define the ensemble-averaged density 



n(r,t) = {n{r,t)), 



(25) 



and current density 



j{r,t)^{j{r,t)), (26) 
where the current operator is defined as 



with 



Pi + eAextiri^t) 



m 



(28) 



the velocity operator of particle i, and the symbol 
{A, B} = {Al3 + 13 A) is the anti-commutator of any two 



The last term on the right-hand side of Eq. ([291) is iden- 
tically zero for bath operators that are local in spacers 
namely 



2V'tn(r, t)V - ytt>n(r, t) - n{r, t)V'<v) = 0. (30) 



Most transport theories satisfy this requirement since 
the action that a true bath does on the system is derived 
from microscopic mechanisms (e.g., inelastic processes) 
which are generally locali^ If this were not the case, 
then this term would represent instantaneous transfer of 
charge between disconnected - and possibly macroscopi- 
cally far away - regions of the system without the need of 
mechanical motion, represented by the first term on the 
right-hand side of Eq. (|29p . This instantaneous "action at 
a distance" is reminiscent of the postulate of wave-packet 
reduction whereby the system may change its state in a 
non-unitary way upon measurement. 

Here, it is the result of the memory-less approximation 
that underlies the stochastic Schrodinger equation (|15p . 
By assuming that the bath correlation times are much 
shorter than the times associated with the dynamics of 
the system (in fact, in the Markov approximation these 
correlation times are assumed zero), we have lost in- 
formation on the microscopic interaction mechanisms at 
time scales on the order of the correlation times of the 
bath. In other words, we have coarse-grained the time 
evolution of our system, and we are therefore unable to 
follow its dynamics on time scales smaller than this time 
resolutioni^ 

In the following, we will assume that the condition (150)1 
is identically satisfied or, if it is not, at any given time, 
given a physical ensemble-averaged current density, a 
unique solution for the ensemble-averaged density can 
be found from Eq. 



E. Equation of motion for the current density 

Similarly, we can derive the equation of motion for the 
ensemble-averaged current densitjii^ 



dtj{r,t) = !^dtA,x,{r,t) - ^ x [V x A.t(r,i)] 



m 



(31) 
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where we have definedsi 

Gir, t) V^j{r, t)V - ij(r, t)V^V - \v^Vj{t, t), 
T[r, t) = - ^5{r - fi)VjU[f, - fj) + mV • a{r, t) 

(32) 

with the stress tensor a (r, t) given by 

^^,j{r,t) = -^^{vtAv3,S{r - fk)}}. (33) 

k 

The first two terms on the rhs of Eq. pTjl describe 
the effect of the apphed electromagnetic field on the dy- 
namics of the many-particle system; the third is due 
to particle-particle interactions while the last one is the 
"force" density exerted by the bath on the system. This 
last term is responsible for the momentum transfer be- 
tween the quantum-mechanical system and the environ- 
ment. 



in an interacting system by a given vector potential, 
then it exists a non-interacting system (the KS system) 
in which we can obtain the same current density by 
applying another suitable vector potential, we will call 
from now on A^f / . 

This is opposed to the general result that an interact- 
ing 1^— representable current density (namely one that is 
generated by a scalar potential V) is not necessarily non- 
interacting representablei^ In particular, it has been 
shown that the mapping between the current density and 
the scalar potential is not invertiblei^i This result shows 
that time-dependent DFT does not necessarily provide 
the exact current density, even if the exact exchange- 
correlation potential is known (albeit it provides the ex- 
act total current for a finite and closed system^^) . With 
some hindsight this is not surprising since there is clearly 
no one-to-one correspondence between a scalar and a vec- 
tor. 

A. The stochastic Kohn-Sham equations 



III. STOCHASTIC TIME-DEPENDENT 
CURRENT-DENSITY FUNCTIONAL THEORY 

Having discussed the physical and mathematical re- 
quirements for the problem we are interested in, we 
can now state the following theorem of stochastic time- 
dependent current-DFT.^" 

Theorem: Consider a many-particle system described 
by the dynamics in Eq. ((2]) with the many-body Hamil- 
tonian given by Eq. ([1]). Let n{r,t) and j{r,t) be 
the ensemble-averaged single-particle density and current 
density, respectively, with dynamics determined by the 
external vector potential Ag^t (?", t) and bath operators 
{Vq}. Under reasonable physical assumptions, given an 
initial condition |^'o), and the bath operators {Vq.}, an- 
other external potential A'^^^ (r, t) which gives the same 
ensemble-averaged current density, must necessarily co- 
incide, up to a gauge transformation, with Aexti^jt). 

The details of the proof of this theorem can be found in 
Ref. Iiol Here, we just mention that the initial condition 
need not be a pure state for the theorem to be valid, but 
may include also the case of mixed initial states. 

The general idea of the proof, following similar ones 
proposed by van Leeuwen^i, and Vignale?^, is to show 
that the external potential A'^^^ is completely deter- 
mined, via a power-series expansion in time, by A^xt, the 
ensemble-averaged current density, the initial condition, 
and the bath operators. 

A lemma of the theorem states that any 
ensemble-averaged current density that is inter- 
acting A— representable is also non-interacting 
A— representable. (A current density is A— representable 
if and only if it can be generated by the application 
of an external potential A.) This implies that if an 
ensemble-averaged current density can be generated 



Let us now assume that we know exactly the vector 
potential A^f / that generates the exact current density in 
the non-interacting system. By construction, the system 
follows the dynamics induced by the SSE (for a single 
bath operator) 

dl'fKs) = [-iHks - \y^y^ \'^Ks)dt -f V\^Ks)dW 

(34) 

where \^ks) is a Slater determinant of single-particle 
wave- functions and 



TT [p^ + eAeff(r,,t)f 



(35) 



is the Hamiltonian of non-interacting particles. 

Note that for a general bath operator acting on many- 
body wave- functions one cannot reduce Eq. ([M]) to a 
set of independent single-particle equations. The rea- 
son is that our theorem guarantees that one can decou- 
ple the quantum correlations due to the direct interac- 
tion among particles, but one cannot generally decouple 
the statistical correlations induced by the presence of the 
environment. These affect the population of the single- 
particle states of the quantum-mechanical system, while 
the quantum correlations are taken into account to all 
orders by the external potential A^ff acting on the KS 
system. It is only when the bath operators act on single- 
particles or on the density that one can write Eq. (j34p 
as a set of equations of motion, one for every KS single- 
particle statcjii 



B. Initial conditions 

The initial condition for the time evolution of the KS 
system has to be chosen such that the ensemble-averaged 
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particle and current densities coincide with those of the 
many-body interacting system. Again, it is important 
to stress that in going from the interacting system to its 
non-interacting doppelganger, the bath operator is not 
modified. On the other hand, the bath operator V gen- 
erally induces transitions between many-body states of 
the interacting Hamiltonian ([T|). Therefore, when rep- 
resented in the non-interacting basis of the KS Hamil- 
tonian it may connect many different single-particle KS 
states. It has been argued that this way the KS system 
will never reach a stationary state even if the coupling 
with the environment is purely dissipativCfH It would be 
thus tempting to modify the bath operator to force the 
KS system into an equilibrium with the external environ- 
mentM- This procedure, however, breaks the theorem we 
have proved, and contains approximations of unknown 
physical meaning. 

In reality, if the true many-body system reaches equi- 
librium with the environment, then the ensemble- average 
current and particle density would attain a stationary 
limit. Since, these are the only two physical quantities 
that the KS system needs to reproduce, the question of 
whether the latter is in equilibrium with the environment 
or not has no physical relevance. 



C. The exchange-correlation vector potential 

The vector potential A^f / (r, t) , acting on the KS sys- 
tem is generally written as the sum of two contributions 



Aeffir, t) = A^^t{r, t) + Ah-xc{r, t), 



(36) 



where A^xt {f, t) is the vector potential applied to the true 
many-body system, and Ah-xc{f, t) is the vector poten- 
tial whose scope is to mimic the correct dynamics of the 
ensemble-averaged current density. From the theorem 
we have proven, Ah-xc{r, t) is a functional of the average 
current density j(r,t'), for t' < t (namely, it is history- 
dependent), the initial condition and the bath operator 
VM 

A common expression would isolate from Ah-xcifji) 
the Hartree interaction contribution from the "rest" due 
to the particle exchange and correlation, namely one 
makes the ansatz 



Ah-xc{r, t) = Ah{r, t) + Axc{r, t) 



where Ah (r, t) is the Hartree contribution to the vector 
potential {t^ is the initial time) 



A,(r,t)= fdt'V [ dr'^^i^. 

Jto J \r-r'\ 



(38) 



The other contribution, Axc{r, t) is again a functional of 
the average current density j{r,t'), for t' < t, the initial 
condition, and the bath operator Vr^ 



In the present case, however, particular care needs 
to be applied to the above ansatz. We have writ- 
ten the Hartree contribution in terms of the ensemble- 
averaged density. This choice, however, requires that the 
exchange-correlation vector potential included also the 
statistical correlations of the direct Coulomb interaction 
at different points in space. These correlations may be 
very large, and possibly much larger than the Coulomb 
interaction between the average densities. The ambigu- 
ity here, compared to the pure-state case, is because in 
a mixed state, quite generally the ensemble-average of 
the direct Coulomb interaction energy contains statisti- 
cal correlations between densities at different points in 
space, namely 



dr / dr' 



,(n(r)) (n(r')) 



J \r-r'\ 

■ (40) 

In actual calculations, one would instead use the form 
of the Hartree potential in terms of the density per el- 
ement of the ensemble. This choice (l38l) makes the KS 
Hamiltonian (j35p stochastic, and therefore, as discussed 
in Sec. [ill no closed equation of motion for the many- 
particle KS density matrix can be obtained. 

Finally, in view of the fact that one can derive a Marko- 
vian dynamics only on the basis of a weak interaction 
with the environment,—!^ as a first approximation, one 
may neglect the dependence of the exchange-correlation 
vector potential on the bath operator, and use the stan- 
dard functional of TD-DFT and TD-CDFT.^-^ [Li]^g ^-^e 
Hartree term, these functionals would also contribute to 
the stochasticity of the KS Hamiltonian ([55)1 .] This seems 
quite reasonable, but only comparison with experiments 
and the analysis of specific cases can eventually support 
it. We thus believe that more work in this direction will 
be necessary. 



IV. CONNECTION WITH A 
DENSITY-MATRIX APPROACH 



From the KS Slater determinants \'^ks)^ solutions 
of the KS equation (|34|) occurring with weight Pi{t) 
iJ2iPiW = 1) in the ensemble, we can construct the 
many-particle KS density matrix [from Eq. (|20p ] 



(37) pKsit) - l*k5W)(*K5WI ^ E ^''W \'^Ksit)){^'Ksit)\- 



Axc{r,t)=Axc j{r,t'),\^o),V 



(39) 



(41) 

This density matrix is, by construction, always positive- 
definite, despite the fact that the KS Hamiltonian and 
possibly the V operator are time dependent. Since in 
general the bath operator acts on many-particle states, 
this many-particle KS density matrix cannot be reduced 
to a set of single-particle density matrices (see Sec.lVlfor 
a numerical ansatz suggested in Ref. to simplify the 
calculations). 

We note first that, in principle, if we knew the ex- 
act functional A^-xc as a functional of the averaged cur- 
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rent density, then the KS Hamiltonian ((35|) would not be 
stochastic and we could derive the equation of motion of 
the many-particle KS density matrix (j4T|) . This equation 
of motion would be Eq. (^5)1 with H replaced by Hks- 

It is important to point out, however, that it is only 
when we start from the stochastic KS equation ([M]) to 
construct the density matrix (|4T|) that we are guaranteed 
that the solution of Eq. ([23]) for the KS density matrix 
maintains positivity at any time. The reverse is not nec- 
essarily true: the equation of motion (p3|) for the KS 
density matrix may, for an arbitrary bath operator V 
or initial conditions, provide non-physical solutions. In 
other words, Eq. (|23p admits more solutions than phys- 
ically allowed, while the SSE always provides a physical 
state of the system dynamics. 

We also stress once more that any approximation to 
Ah-xc, makes the KS Hamiltonian stochastic, namely one 
Hamiltonian for each element of the ensemble, thus mak- 
ing the density-matrix formalism of limited value. In fact, 
by insisting on using Eq. (|23p with these approximations 
would amount to introducing uncontrollable approxima- 
tions in the system dynamics which entail neglecting im- 
portant statistical correlations induced by the bath (see 
also discussion in Sec. IVI|1. 

To see this point explicitly, let us consider the equa- 
tion of motion for an arbitrary operator acting on the 
KS system that evolves according to Eq. (fT7|) with Hks 
replacing H 



d{A) 



i[A, Hks] - ^ [v^VA + AV'^V - 2VAV'^^ \ dt 



{V''A + Av)dW. 



(42) 



Now we take the ensemble average of this equation in 
order to obtain the equation of motion for the ensemble- 
averaged quantities. However, since now Hks is a 
stochastic Hamiltonian, then the ensemble average and 
the commutator between A and Hks do not commute, 
i.e., in general we expect that 



[A, Hks] ^ [A, Hks], 



(43) 



which implies that 



dtPKS ^ - i Pks, Hks 



\ {y^VpKS + PKSV^V ~ 2VpKsV^) , 



(44) 



namely there is no closed equation of motion for the KS 
density matrix. 

In fact, the correct procedure is to evolve the system 
for every realization of Hamiltonians, and then average 
over these realizations, which is what a solution of the 
SSE dMD would provide. 



V. NUMERICAL SOLUTION OF THE SSE 
A. Finite difference equation 

We now discuss practical implementations of the 
SSE ([M]) . First of all we realize that, in going from the 
differential equation (p4)) to a finite difference equation 
that can be solved on a computer, one has to bear in mind 
that dW is on the order of y/dt. Then one has to expand 
the equation for the finite differences and keep terms of 
second order in dW that correspond to first order in dt. 

Here we write down the correct finite difference equa- 
tion starting from In the following we assume that 
the state vector ks is a regular function of time, posi- 
tion, and the Wiener process W , i.e., we assume that the 
derivatives 



5* 



KS 



KS 



dW ' dWdW 



(45) 



exist and are regular. 

Let us define dt — t — t' a small time interval over 
which we integrate the equation of motion for 'i! ks- If we 
expand in series the increment d^KS = ^ Ks{i) — '^ Ks{t') 
we have, 



KS 



— dW + -- — - — dWdW ' 

ui dW 2 dWdW 



{^"dt ' 2 



A direct term-by-term comparison with Eq. (|34p tells us 
that there is a correspondence between d/dW and V so 
that a finite difference scheme can now be implemented 
such that the equation of motion 



KS — 



KS 



\ dt 



At 



KS 



V dW 



AW 



1 



-IHks - ^V^V - *KsAt 



+ V^KS^W, 



(47) 



is the correct first-order equation in Ai. Eq. (|T7)) can now 
be solved by a variety of different numerical techniques, 
and we refer the reader to other work for a discussion of 
such methods. ^i^i^i^^ 

The important point is that one evolves these equa- 
tions in time for every realization of the stochastic pro- 
cess and then averages over the different realizations (in 
Sec. IVII we give an explicit example of such calculation 
showing the convergence of the results with the number 
of realizations). 



B. The non-linear SSE 

The norm of the state vector solution of the SSE ([M)) is 
preserved on average but not for every realization of the 
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stochastic process.-^'J^ This may slow down the conver- 
gence of the results as a function of the number of real- 
izations of the stochastic process. It is thus more conve- 
nient to solve a non-linear SSE which gives an equivalent 
solution as the linear SSE p4|) . This can be easily done 
by first calculating the differential [in the Ito sense ([7])] 
of the square modulus of \'^ks) 



= 2R\\^Ks\\^dW, 
where we have defined 



R 



2 \\^KS\ 

By using the power expansion 
dW'^KsW 



(48) 



(49) 



d^W^KsW 
1 



1 



KS\ 



KS 



|2)3 



d\\^Ks\?d\\^Ks\?m 



The finite-difference equation for this case is 



9$ 



KS 



dt 



At 



KS 



V dW 



we can derive the equation of motion for the state vector 
normalized at every realization of the stochastic process: 



1$ 



I* 



KS) 



KS/ 



I* 



KS\ 



(51) 



which is (see also Ref. 02 



KS) — 



RV - -Rn 



-{V - Ri)\<i>Ks)dW. 



\^Ks)dt 
(52) 



This non-linear equation of motion, by construction, is 
equivalent to the linear SSE ([51)1 . 



AW 



= I -iHks - \v^V - + RV- 1-RH ] $KsAt + {V- Ri)<i>KS^W. 



(53) 



C. Single- particle order- A'^ scheme 

Due to the presence of the environment, it is still a 
formidable task to solve the equations of motion of S- 
TDCDFT for arbitrary bath operators. In fact, as we 
have already discussed, the bath operators generally act 
on Slater determinants and not on single-particle states. 
If we have N particles and retain M single-particle states, 
this requires the solution of — 1 elements of the state 
vector (with Cjjf = Ml/Nl{M - Ny. and the -1 comes 
from the normalization condition). In addition, one has 
to average over an amount, call it m, of different realiza- 
tions of the stochastic process.''^ The problem thus scales 
exponentially with the number of particles. 

However, it was recently suggested in Ref. [2^ that 
for operators of the type A = Aj , sum over single- 
particle operators (like, e.g., the density or current den- 
sity), the expectation value of A over a many-particle 
non-interacting state with dissipation can be approxi- 
mated as a sum of single-particle expectation values of 
Aj over an ensemble of N single-particle systems with 
specific single-particle dissipation operators. In particu- 
lar, the agreement between the exact many-body calcula- 



tion and the approximate single-particle scheme has been 
found to be excellent for the current density,^ We refer 
the reader to Ref. for the numerical demonstration of 
this scheme and its analytical justification. The physical 
reason behind it is that, due to the coupling between the 
system and the environment, highly-correlated states are 
unlikely to form. 

Here, for numerical convenience, we will adopt the 
same ansatz which in the present case reads, 



N 



{•^KS\A\^KS) 



«s\A,\^i's 



(54) 



with \<j>''j(s) single-particle KS states solutions of 



^'ks) 



^ {p + eA,ff{r,t)f ly.y 

+Vsp\^'Ks)dWit), 



\<l>iis)dt 
(55) 



with Vsv an operator acting on single particle states (see 
Refs. [2a and[ll| and next section for explicit examples of 
such operator) 4^ 
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For convenience, also in the present case we can nor- 
malize the single-particle KS states for every realization 
of the stochastic process by defining 



KS/ 



KS' 



i^'Ksl 



(56) 



and thus solve the non-linear SSE 



d\4>- 



KS' 



. (p + eAe//(r, t)f 
' 2m 



KS 



1 

sp" sp 



where 



HVsp^Rjm'Ks)dWit), 

li^^siivJp+yspM^s) 



110- 



(57) 



(58) 



KS\ 



The discretization of these equations is then done sim- 
ilarly to what we have explained in the previous section. 



account correctly the statistical correlations induced by 
the bath, while the SSE naturally accounts for the the 
stochasticity introduced in the Hamiltonian by the in- 
teraction potential. In fact, we find that when both the 
initial and final state are pure, both approaches provide 
the same equilibrium state. However, the correspond- 
ing dynamics are different. In particular, the relaxation 
time obtained from the evolution of the density matrix 
is shorter than the relaxation time obtained from the av- 
erage over many realizations of the dynamics obtained 
from the SSE. 

The differences between the two approaches are even 
more striking when we consider the evolution towards 
a state that contains at least two major contributions 
coming from different states. In this case, also the final 
steady states obtained from the density matrix and the 
SSE are different. 

These cases exemplify what we have discussed all 
along: if one insists on using a closed equation of mo- 
tion for the density matrix of the type (|23p with stochas- 
tic Hamiltonians, uncontrolled approximations are intro- 
duced which lead to an incorrect dynamics. 



VI. AN EXAMPLE: A GAS OF LINEAR 
HARMONIC OSCILLATORS 

Stochastic-TDCDFT has been applied to the study of 
decay of exited He atoms and its connection to quantum- 
measurement theory..^ ^ It can describe the dynamics of 
bosons as well. In this section, we apply it to the anal- 
ysis of the dynamics of an interacting ID Bose gas con- 
fined in a harmonic potential and coupled to a uniform 
external environment that forces the gas towards some 
steady state. Since neither the external potential nor 
the bath are time-dependent, we expect that the boson 
gas reaches a steady state configuration when coupled 
with the uniform external bath. Finally, we assume that 
the bath forces the system towards certain eigenstates of 
the instantaneous interacting boson Hamiltonian. The 
bosons are interacting via a two-body contact potential, 
i.e., Uint(x,x') oc 5(x — x'). This potential correctly 
describes the important case of Alkali gases in which 
the Bose-Einstein condensation has been experimentally 
observed i^Si^iH 

The purpose of this section is to compare the dynam- 
ics of the boson gas obtained from the SSE [Eq. p^ ] and 
the quantum master equation [Eq. For this reason 

the value of the physical parameters (the strengths of the 
confining potential, the particle-particle interaction, and 
the system-bath coupling) is arbitrary and chosen only 
for the sake of this comparison. We will report elsewhere 
a more realistic study of the dynamics of this important 
physical system4^ When no interaction between particles 
is included both approaches are clearly equivalent. How- 
ever, when interactions are included, the Hamiltonian of 
the system becomes stochastic and, as previously dis- 
cussed, the quantum master equation does not take into 



A. Macroscopic occupation of the ground state 

We begin with the study of the dynamics of the macro- 
scopic occupation of the ground state induced by energy 
dissipation towards the degrees of freedom of an exter- 
nal bath. The external bath forces the system to reach 
a state of zero temperature or minimal free energy, i.e., 
the ground state of the Hamiltonian. One possible form 
of this bath operator is, in a basis set that makes the 
Hamiltonian diagonal at each instance of time^^i 



V=5 



/Gill 




(59) 



\ ... ) 



where <5 is a coupling constant with dimensions of the 
square root of a frequency [we set 5 = ^/uni in what fol- 
lows, with Wo the frequency of the harmonic confining 
potential - see Eq. (|63p ] . We do not expect that this op- 
erator fulfills Eq. ([50]) since in a real-space representation 
it would allow for the localization of the particles without 
an effective current between two distinct points in space. 
In this section, however, we are more interested in the 
kind of dynamics this operator generates in our quantum 
system, and the comparison with the dynamics obtained 
from the quantum master equation. We expect, indeed, 
that the condition Eq. ([5n|) is violated both in the SSE 
and in the quantum master dynamics. 

The operator (j59p mimics the energy dissipation in the 
system, with the external bath absorbing the bosons' ex- 
cess energy and cooling down the boson gas. One could 
argue that this is the generalization to the many-state 
system of the bath considered in Ref . frj) . We can thus 
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conclude that the effective temperature of the bath we 
consider here is zero. 

The Hamiltonian of the boson system (in second quan- 
tization) when the bath is not present reads 

dxdx' ijjHx)il^\x')Uint{x - x')V'(a;')V'(aOgO) 



coupled to the boson gas, 



where '4){x) destroys a boson at position x, Vext{x) is a 
confining potential, and Uint{x — x') is the boson-boson 
interaction potential. For dilute boson atomic gases the 
interaction potential can be substituted with the contact 
potential, i.e., 

Uint{x - x) = go{N - l)5{x - x') = g5{x - x') (61) 

where is determined by the scattering length of the 
boson-boson collision in the dilute gas, and TV is the total 
number of bosons in the trap, so that ||'0|| = li^ 

With standard techniques, and in the Hartree approx- 
imation, we can go from the equation of motion for the 
annihilation operators to the equation of motion for the 
state of the system ^'(a;, i), when the external bath is not 



idt^{x,t) = 



^{x,t)+gn{x,t)'ii{x,t) 
(62) 

where n{x,t) — \'^(x,t)\'^ is the single-particle density of 
the boson gasi^^i^ 

Equation ([62|) (and its generalization to 2 and 3 
dimensions) has received a lot of attention since it 
correctly describes the dynamics of a Bose-Einstein 
condensate.i^2i2ii2^ 

In the following we will focus on the case of a ID har- 
monic confining potential, i.e.. 



Vext{x) = ^muigx'^. 



(63) 



A harmonic confinement is created, e.g., in the magneto- 
optical traps used in the experimental realization of 
the Bose-Einstein condensation of dilute boson Alkali 
gases J^2i2ii2^ 

When the boson system is coupled to the external envi- 
ronment, we assume that the Hamiltonian is not affected 
by the coupling and the state of the system '^{x, t), that 
is now stochastic, evolves according to the SSE 



J 



d-^ix^t) 



7r~T^ + 7;^^o^'^ +gn{x,t) ) ^{x,t)dt 
2m dx^ 2 



-V^V'i'ix, t)dt + V^{x, t)dW 



(64) 



where this equation of motion has to be interpreted in accordance to the discussion of the previous Sections. For 
numerical convenience, we rescale this equation in terms of the physical quantities Wq, Xq — l/^mwo, and g = g/xo 
to arrive at 



d^'(x,i) = -iujQ { - — + ^ + —n(x,t)xo] ■i>(x,t)dt- -V''V'i>(x,t)dt + V-i'(x,t)dW. 
' 2 dx'' 2xq ujo / 2 



(65) 



r 



We begin by considering the case of non-interacting 
bosons, i.e., we set g — 0. In this case, the Hamilto- 
nian admits a natural complete basis, the set of Hermite- 
Gauss wave-functions 



ipj{x) = 



zHj{x/xo)e 



(66) 



where the polynomials {Hj} satisfy the recursion relation 



Hj+i(x) — 2xHj(x) — Hj^i 



(67) 



and Ho{x) = 1, Hi{x) = 2x. If we expand the wave- 
function "^[x,t) = cij{t)ipj{x), and make use of the 



orthonormality properties of the Hermite-Gauss wave- 
functions, we obtain the (stochastic) dynamical equation 
for the coefficients a, , 



(68) 



dttj 



where Hij — (j + l/2)ujQdij and V is given by Eq. ((59| . 

Together with Eq. (|68p we can study the dynamics 
of the density matrix via the quantum master equation 
Eq. which in the same spatial representation as 

Eq. jMl) reads 



ytPij 



= E i^ikPk] - PikHkj) + E ( ^'ikPkk'Vl 



k,k' 



^V^k^kk'Pk'j 



^PikVkk'Vk'j 



(69) 
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The connection between Eq. (ISS)) and Eq. is estab- 
lished by the identity pij = a*aj vahd for any pair of 
indexes i and j. We solve Eq. ([68|) numerically with a 
4th order Runge-Kutta evolution scheme, after we have 
mapped the dynamics to its norm-preserving equivalent 
form (see Sec.lVland the discussion therein) i^ii^ For con- 
sistency we solve the master equation with a 2nd order 
Runge-Kutta evolution scheme (in fact with the more 
refined Heun's scheme) <^ 

We report in Fig. ^ the dynamics of the probability 
of occupation of the ground state, po — |aoP — poo for 
various realizations of the stochastic field in Eq. to- 
gether with the dynamics obtained from the evolution of 
the density matrix l|69p. Here, we have included the first 
20 levels of the free Hamiltonian, and we have chosen as 
initial condition a2o(0) = 1 and set the other coefficients 
to zero. We have set the mass of the particles to 1, and 




FIG. 2: (Color online) Non-interacting bosons - Occupation 
probability of the ground state versus time calculated from 
the evolution of the state via the SSE averaged over differ- 
ent runs (1, 10, 50 and 100) and via the equation of motion 
for the density matrix for the non-interacting boson case. It 
is evident that with only 50 realizations the agreement be- 
tween the SSE and the density matrix equation is excellent. 
In the inset we show the relative difference between the two 
dynamics for 100 realizations of the stochastic process. 

used a time step ujoAt = 20/2^^ = 6 x 10"'*. A further 
decrease of this time step does not affect the results sig- 
nificantly. From Fig. [21 it is evident that when we collect 
a large enough statistics the results of the SSE and the 
master equation coincide for the non-interacting boson 
case: Already for 50 runs of the SSE the difference be- 
tween the two dynamics almost vanishes.^ In the inset 
of Fig. [2] we report the relative difference between the 
occupation number of the ground state with the two dy- 
namics, Ipo™ — Po'^'^l/Po'*^- We see that this difference, for 
100 runs, is generally lower than 5%, a quite satisfactory 
result. 

In Fig.[3l we report the density profile for the system at 
different instances of time obtained from the SSE. Start- 
ing from a pure state, where the highest energy state 



is occupied (panel c), the system relaxes towards the 
ground state. As it is clear from panel a) in Fig. [3] the 
system, at tujQ = 20, still occupies certain high energy 
states [see, e.g., the tail at x < of panel a)]. 
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FIG. 3: Non-mteracting bosons - Plot of the averaged den- 
sity profile, n{x) x Xq, for various instances of time calculated 
from the SSE. The system evolves from the maximum occu- 
pation of the highest excited state (panel c) , to the maximum 
occupation of the ground state. We have averaged over 100 
realizations of the stochastic process. 

We now turn on the particle-particle interaction Uint- 
This corresponds to adding to the free Hamiltonian Hi^j 
an interaction part, iJ™* that in the basis representation 
of the Gauss-Hermite polynomials reads 

Hr:f ^Y.^^,J,k,,ala, (70) 

k,q 

where Fij.k,q is the 4th-rank tensor defined as 

/oo 
dx H,{x)Hj{x)Hk{x)Hg{x)e-^'=\ (71) 
-oo 

A long but straightforward calculation brings us to an ex- 
plicit expression of Fij-k,q in terms of Euler gamma func- 
tions and a hypergeometric functioui^i^i It can be shown 
that the hypergeometric function reduces to the summa- 
tion of a few - at most min(z,j) - terms. In the case of 
the density matrix approach the interaction Hamiltonian 
is immediately written as 

^"* = E^.^;'=-'?'°M- (72) 

k,q 

In solving the dynamics of the system described either 
by the SSE (|65p or the master equation ((69)l . we have 
assumed that at any instance of time the bath operator 
brings the system towards the instantaneous ground state 
of the interacting Hamiltonian Hi j + H™-*- . 

In Fig. [4] we plot the occupation probability pj [t) of 
the state j for the first 3 levels of the free Hamiltonian 
{pj{t) = |aj(t)p from the SSE or Pj(t) = Pjj{t) from 
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the density matrix). We have assumed an interaction 
of strength g/ujQ = 5, and a time step WpAt = 60/2^'' 
and we have performed 100 independent runs of the SSE. 
While it evident that the system reaches the same steady 
state according to the two equations,'*^ it is also clear 
that the state calculated with the SSE relaxes slower 
than the state obtained from the density matrix equation. 
This is a spurious effect in the density matrix dynamics 
where the average density defines the interaction poten- 
tial. This does not take into account the fluctuations of 
the state, and hence of the stochastic Hamiltonian. 

We have also tested that the steady state reached dur- 
ing the dynamics is consistent with the theory of the 
eigenstates of the Gross-Pitaevskii equation4Si22, In par- 
ticular, the ground state of the interacting system, when 
the interaction is strong, can be obtained by neglecting 
the kinetic contribution to the Hamiltonian. In this case, 
a good approximation to the ground state density reads 

9-^0 

-l-terms proportional to (73) 

where /i, the chemical potential, is determined by the nor- 
malization condition, and 9(x) = if x < and 9{x) = 1 
if cc > 0. 

In the inset of Fig. |4] we plot the density obtained at 
tu!o — 60 from the SSE (black, dashed line) and the den- 
sity obtained from the approximation (|73p (orange, solid 
line). Notice that the value of the parameters g and /i 
have been obtained from the best fit with the numerics: 
indeed one can show that the approximation (|73p is ex- 
act in the limit of very large interaction ;^^'^^ which is not 
reached in our calculations. 

In Fig. [5] we report the value of the ground state energy 
of the interacting Hamiltonian versus time as calculated 
from the SSE and the master equation. Again the dif- 
ference between the relaxation times calculated from the 
two dynamics is evident. In the inset of Fig. [5] we report 
the energy of the first excited state. 

To summarize this section, we have described the dy- 
namics of the relaxation of a confined ID boson system 
towards the ground state induced by a given external 
bath. The final state we have obtained is consistent with 
the eigenstate of the ID Gross-Pitaevskii equation. Our 
main result is that, although the SSE and the master 
equation reach the same final state, the dynamics de- 
scribed by these equations show important differences, 
and physical quantities, like, e.g., the relaxation time, 
differ. In particular, the density matrix approach, which 
at any instant of time employs the average density to 
construct the interaction Hamiltonian, underestimates 
the fluctuations induced by the bath on the stochastic 
Hamiltonian. These fluctuations are correctly taken into 
account in the SSE. 




FIG. 4: (Color online) Interacting bosons - Occupation prob- 
ability of the first 3 lowest energy levels of the non-interacting 
Hamiltonian versus time calculated via the SSE (black, solid 
lines) averaged over 100 independent runs, and via the equa- 
tion of motion for the density matrix (red, dashed lines) . The 
time it takes the system to reach steady state is different for 
the density matrix approach and the SSE, with the former 
underestimating the relaxation time. This is due to the in- 
clusion in the master equation of the average density in the 
interaction potential, thus neglecting important fluctuations 
that can slow down the relaxation dynamics. In the inset 
we compare the equilibrium density (black, dashed line) with 
the one obtained from the Thomas-Fermi approximation to 
the ground state (orange, solid line). 
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FIG. 5: (Color online) Interacting bosons - Time evolution of 
the energy of the ground state of the Gross-Pitaevskii Hamil- 
tonian as calculated from the SSE (black, solid line) and mas- 
ter equation (red, dashed line) . The final value of the energy 
is the same, but the relaxation dynamics is different in the 
two formalisms with the master equation considerably under- 
estimating the relaxation time. In the inset we report the 
dynamics of the first excited state of the Gross-Pitaevskii 
Hamiltonian obtained from the SSE (black, solid line) and 
the master equation (red, dashed line). 



B. Competition between states 

Let us now consider the more common case in which 
the environment drives the system toward a mixed steady 
state. To simplify the discussion we consider only three 
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FIG. 6: (Color online) Non-interacting bosons - The occupa- 
tion probabilities for a three level system calculated via the 
SSE (|65[) for non- interacting bosons {g = 0). The results 
obtained via the master equation for the density matrix are 
indistinguishable on this scale from those obtained with the 
SSE. In the inset we show the difference between the occu- 
pation probabilities of the lowest energy level calculated from 
the SSE and the master equation. This difference is in mod- 
ulus lower than 5 x 10"'^ at any instant of time. 



single-particle levels and the bath operator forces the sys- 
tem towards two different states. We choose, in a basis 
in which the Hamiltonian is diagonal, the operator 



V 




(74) 



i.e., the operator drives the system, with equal strength, 
towards the lowest and highest energy levels of the inter- 
acting Hamiltonian. As we will see, the final state is a 
superposition of these two states with a significant contri- 
bution coming from the middle level. At first glance this 
might seem surprising. However, we have to remember 
that, e.g., in the quantum master equation, the equilib- 
rium states are determined by the kernel of the super- 
operator. This super-operator contains powers of the op- 
erator V, that in turn contains a finite contribution from 
the middle level. A similar reasoning applies to the SSE. 

To begin with our analysis of this system, we consider 
the non-interacting case g = 0; we set as before S = ^/uJo^, 
and we start from the fully occupied highest energy level, 
i.e., 03(0) = 1. In Fig. [5] we plot the occupation prob- 
abilities for the 3 levels calculated via the SSE (p5|) . In 
this case, to reduce the stochastic noise even further, we 
have performed 1000 independent runs of the SSE and 
used, in both dynamics, usoAt — 20/2^*. As we can see 
from Fig. O at steady state the bath operator forces the 
system to occupy the lowest and the highest energy levels 
with equal probability, while a finite occupation proba- 
bility of the middle level appears. This mixing prevents 
the system to reach a pure steady state and some finite 
correlation between the energy levels, that appears for 
example in the finite off-diagonal elements of the density 
matrix, persists in the long-time regime. 

Again, for this non-interacting case the dynamics ob- 
tained from the SSE and master equation are indistin- 
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FIG. 7: (Color online) Interacting bosons - Plot of the dy- 
namics of the occupation numbers of the lowest and high- 
est energy level, po and P2, respectively, calculated from the 
SSE (|65[) (black, solid line) and the master equation for the 
density matrix (|23p (red, dashed line). Panel a), g/uJo = 0.1: 
for small interaction the two dynamics show small differences. 
Panel b), guio = 0.5: for intermediate interaction strength the 
differences between the two dynamics is a large fraction of the 
occupation number. Panel c), for large interaction g/ujo ~ 1 
the difference between the two dynamics for the lower level 
decreases. In this particular case, this is due to the presence 
of the second energy level (not shown in the figure) that is 
little affected by the interaction. 



guishable on the scale of the plot of Fig. [6j^ In the inset 
of Fig. we report the difference between the ground 
state occupation probability as calculated from the SSE 
and from the density matrix approach. This difference 
is, in amplitude, smaller than 5 x 10^'^, and by increas- 
ing the number of independent runs, it decreases. To test 
our numerical code, we have also compared the numerical 
solution with the exact dynamics obtained from the ana- 
lytical solution of the master equation (which is feasible 
because we have only three states). Since the numerical 
and analytical solutions are essentially the same, we do 
not find necessary to report the analytical solution here. 



We now turn on the particle-particle interaction (|6ip . 
Fig. [7] reports the time evolution of the occupation 
number of the lowest and highest energy levels of the 
free Hamiltonian for different strengths of the particle- 
particle interaction. As expected, the interaction opens 
a gap in the occupation numbers between the highest and 
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lowest energy levels. Most importantly, we see that for 
intermediate values of the interaction the steady states 
calculated with the SSE and the master equation differ. 
This difference is not monotonic with the interaction, and 
it is state-dependent. We see indeed that for relatively 
strong interaction g/ioa > 1, this difference is smaller 
than for g/ujQ — 0.5; more so for the lowest state than 
the highest one. This is due to the fact that the middle 
energy level (not shown in the figure), which is almost 
unaffected by the variation of the interaction strength 
and whose dynamics is almost the same for the SSE and 
the master equation, "blocks" the transformation of the 
highest energy level to low occupation numbers. For very 
strong interaction g/ujQ = 5 (not shown in the figure), the 
occupation numbers calculates via the SSE and the den- 
sity matrix approach, almost coincide. 

The above example shows that when the bath drives 
the system towards a mixed state, also the final states 
(not just the dynamics) obtained from the density ma- 
trix according to the master equation (^5)1 and the SSE 
may be different. In the particular case considered here, 
this is due to the fact that the final state is sensitive to 
the frequency of the confining potential (as can be shown 
with the exact analytical solution of the non-interacting 
system). The SSE and the master equation create differ- 
ent effective interaction potentials that renormalize the 
frequency of the confining harmonic potential. This dif- 
ferent renormalization shows up in the different steady 
states. This important difference is again due to the fact 
that in the master equation the interactions are included 
using the average particle density, thus neglecting the 
true stochasticity of the Hamiltonian. Small differences 
in the effective potential (confining plus interaction) thus 
results in macroscopic differences in the steady states. 
The fact that the dynamics of the interacting system de- 
scribed by the master equation ([^5]) is so sensitive to the 
interaction potential and does not reproduce correctly 
the dynamics and/or the steady states of the system un- 
dermines the applicability of an equation of motion for 
the density matrix to the stochastic extension of TDDFT 
and TDCDFT. 



VII. CONCLUSIONS 

In this paper, we have discussed in detail a functional 
theory of open quantum systems we have named stochas- 
tic TDCDFT. This theory, based on a theorem we have 
previously proved in Ref. [13, extends DFT to the dy- 
namical interaction of quantum systems open to exter- 
nal environments, when the latter satisfy a memory-less 
dynamics. The starting point of the theory is a stochas- 
tic Schrodinger equation for the iV— particle state vec- 
tor, which provides a conceptually transparent way of 
describing open quantum systems. 

We have discussed the mathematical assumptions of 
the theory, the numerical solution of the corresponding 
equations of motion, and compared it to a possible for- 



mulation in terms of a density-matrix approach based on 
quantum master equations. We have shown that due to 
the dependence of the KS Hamiltonian on microscopic 
degrees of freedom, and its time-dependence, a density- 
matrix approach to a stochastic DFT is not a solid alter- 
native to this problem. In fact, due to these conditions, 
there is not necessarily a closed equation of motion for 
the density matrix, and if one insists on using a quantum 
master equation, the solutions of such an equation may 
not be physical for all cases. 

As an example of application, we have used the the- 
ory to study the dynamics of a ID gas of excited bosons 
confined in a harmonic potential and in contact with an 
external bath. This is a problem previously inaccessi- 
ble by standard DFT. Along similar lines, we expect this 
theory to find application in a wide range of problems 
where DFT methods could not be applied, such as en- 
ergy transport and dissipation, dephasing induced by an 
environment, quantum measurement and quantum infor- 
mation theory, phase transitions driven by dissipativc ef- 
fects, etc. 

From here, an interesting (and non-trivial) exten- 
sion of stochastic TDCDFT would be to environments 
with finite auto-correlation times. This leads to non- 
Markovian dynamics with memory kernels and more 
complex stochastic Schrodinger equations liiiii If a sim- 
ilar theorem as that we have demonstrated here can be 
proved for these cases as well, we could study an even 
larger class of open quantum system problems, where 
memory effects in the bath dynamics are of particular 
importance. 

Another possible extension of the theory would be to 
investigate the noise properties of the quantum system. 
This would provide even more information on the system 
dynamics. An extension of S-TDCDFT to this problem 
seems possible but not trivial. The reason is because the 
noise is an n-time correlation function (where n indicates 
the moments of the observable), and as such it cannot 
be written simply in terms of the expectation value of 
an observable. It is thus not obvious what is the phys- 
ical variable conjugated to the noise of given moment. 
One could clearly calculate the moments of the current 
using the present form of S-TDCDFT. How good this 
approximation is compared to the exact noise (even if 
one knows the exact functional of S-TDCDFT) is an is- 
sue that, like other applications of DFT beyond its basic 
theorems (e.g., the assignment of a physical meaning to 
the KS states) , must be addressed at an "empirical" level 
by comparing with experiments or available analytical re- 
sults. 

Finally, another important direction of study would be 
the development of functionals in the presence of baths. 
Clearly, this cannot be done for arbitrary baths, and spe- 
cific cases, such as a bath of harmonic oscillators, would 
be a good starting point. It would be interesting to know 
if an approximate functional with a clear physical inter- 
pretation can be obtained, and how different it is from 
the functionals in the absence of bath interaction. Un- 
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til then, the best we can do is to apply the available 
functionals, justify their use on the basis of the weak in- 
teraction between the system and the environment, and 
compare the results with available experimental data or 
analytical results. 
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This assumption is needed in the proof of the theorem of 
Stochastic- TDCDFT, see Ref [TO. 

In the following, see Sec.|Vl we will derive the finite difi'er- 
ence equation that is satisfied by the wavefunction ^ . 
*i In Eq. (|32|) . Vj contains the derivatives with respect to 
the coordinates of the j-th particle, i.e., in 3D Vj = 

{d:^,,dy^,d^.). 

*2 A density-matrix formalism would be even more computa- 
tionally demanding, requiring the solution of {C^ -I- 2) x 
{C^ — l)/2 coupled differential equations, even after taking 
into account the constraints of hermiticity and unit trace 
of the density matrix. 

Note that the theorem of S-TDCDFT is still valid, and 
Eq. (|54|l would be exact (and not an approximation) , if we 
choose the bath operators to act on single-particle states 
(or the density) to begin with. 
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In calculating the time evolution with the SSE we make 
use of the techniques discussed in Sec. [V] 
We expect that, for non-interacting particles, the devia- 
tion between the dynamics obtained via the density matrix 
equation and the SSE, scales as if m is the number 

of independent runs on which we average the SSE. 
The initial state is pure and the bath is selecting only a par- 
ticular state thus forcing the system towards another pure 



state. Moreover, we can prove that if the system evolves 
from the ground state, the stochastic part vanishes on this 
state, and then the boson gas remains in the ground state 
of the interacting Hamiltonian. 

Only the dynamics obtained from the SSE is reported in 
Fig.H 



